Phase Transition in Heisenberg Stacked Triangular Antiferromagnets: End of a 
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By using the Wang-Landau flat-histogram Monte Carlo (MC) method for very large lattice sizes 
never simulated before, we show that the phase transition in the frustrated Heisenberg stacked 
triangular antiferromagnet is of first-order, contrary to results of earlier MC simulations using old- 
fashioned methods. Our result lends support to the conclusion of a nonperturbative renormalization 
group performed on an effective Hamiltonian. It puts an end to a 20-year long controversial issue. 
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I. INTRODUCTION 

When a spin cannot fully satisfy energetically all the 
interactions with its neighbors, it is "frustrated". This 
situation occurs when the interactions are in competi- 
tion with each other or when the lattice geometry does 
not allow to satisfy all interaction bonds simultaneously 
as seen for example in the triangular lattice with an anti- 
ferromagnetic interaction between the nearest-neighbors. 
Effects of the frustration in spin systems have been ex- 
tensively investigated during the last 30 years. Frus- 
trated spin systems are shown to have unusual properties 
such as large ground state (GS) degeneracy, interesting 
GS symmetries, successive phase transitions with com- 
plicated nature, partially disordered phase, reentrance 
and disorder lines. Frustrated systems still constitute 
at present a challenge for theoretical, experimental and 
simulational methods. For recent reviews, the reader is 
referred to Ref. [l|. 

The nature of the phase transition in strongly frus- 
trated spin systems has been a subject of intensive in- 
vestigations in the last 20 years. Theoretically, these 
systems are excellent testing grounds for theories and 
approximations. Many well-established methods such as 
renormalization group (RG), high- and low-temperature 
series expansions etc often failed to deal with these sys- 
tems. Experimentally, data on different frustrated sys- 
tems show a variety of possibilities: first-order or second- 
order transitions with unknown critical exponents etc. 
(see reviews in Ref. |l|). One of the most studied sys- 
tems is the stacked triangular antiferromagnet (STA): 
the antiferromagnetic (AF) interaction between nearest- 
neighbor (NN) spins on the triangular lattice causes a 
very strong frustration. It is impossible 1 - to fully sat- 
isfy the three AF bond interactions on each equilateral 
triangle. The GS configuration of both Heisenberg and 
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XY models is the well-known 120-degree structure. The 
cases of XY (N = 2) and Heisenberg (N = 3) spins on 
the STA have been intensively studied since 1987. For de- 
tails, see for example the review by Delamotte et al£. Let 
us briefly recall here some main historical developments. 
Kawamura^ has conjectured by a two-loop RG analy- 
sis and Monte Carlo (MC) simulations that the transi- 
tion in XY and Heisenberg models belong each to a new 
universality class in dimension d = 3. Since then there 
have been many other calculations and simulations with 
contradictory results. For example, Azaria et a£ sug- 
gested from a non-linear sigma model that if the transi- 
tion is not of first order or mean-field tricritical then it 
should be 0(4) universality. Numerical simulations^^ 
however did not confirm these conjectures. Antonenko 
et al£ went further in a four-loop RG calculation with 
a Borel resummation technique. They concluded that 
the transition is of first order. From 2000, Tissier and 
coworkers 10 i 1:L i 12 have carried out a nonperturbative RG 
study of frustrated magnets for any dimension between 
two and four. They recovered all known perturbative 
one-loop results in two and four dimensions as well as 
for the infinite spin-component number N — > oo. They 
determined N c (d) for all d and found N c (d = 3) = 5.1 
below which the transition is of first order in contradic- 
tion with the conjecture of the existence of a new chiral 
universality class by Kawamurai^ They explained why 
theories and simulations have encountered so far many 
difficulties by the existence of a whole region in the flow 
diagram in which the flow is slow: the first-order char- 
acter for N = 2, 3 is so weak that the transition has 
a second-order aspect with "pseudo" critical exponents. 
They calculated these pseudo exponents and found that 
they coincided with some experimental data. While this 
scenario is very coherent, we note that in this nonpertur- 
bative RG technique, the real Hamiltonian is truncated 
at the beginning and replaced by an effective one. How- 
ever, as will be seen in this paper, the nonperturbative 
results are well confirmed. 

Let us recall some results on the XY case. Early MC 
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results on XY STA have been reviewed by Loison^ Un- 
til 2003, all numerical simulations found a second-order 
transition with exponents. A numerical breakthrough 
has been realized with the results of ItakuraM who used 
an improved MC renormalization-group scheme to in- 
vestigate the renormalization group flow of the effective 
Hamiltonian used in field-theoretical studies for the XY 
STA. He found that the XY STA exhibits a clear first- 
order behavior and there are no chiral fixed points of 
renormalization-group flow for N=2. In 2004, Peles et 
ali£ have used a continuous model to study the XY STA 
by MC simulation. They found evidence of a first-order 
transition. In 2006, Kanki et ali£, using a microcanoni- 
cal MC method, have found a first-order signature of the 
XY STA. While these recent simulations have demon- 
strated evidence of first-order transition for the XY STA 
in agreement with the nonperturbative RG analysis, all 
of them suffered one or two uncertain aspects: the work 
of Itakura has used a truncated Hamiltonian, the work 
of Peles et al has used standard MC methods and the 
work of Kanki et al used a traditional microcanonical MC 
technique. Using a very high-performance technique for 
weak first-order transitions, the so-called Wang-Landau 
flat-histogram method^ we have recently carried out 
simulations on the XY STA. We have found clearly a 
first-order transition in that system confirming results 
of other authors and putting an end to the controversy 
which has been lasting for 20 years. 

For the Heisenberg case, ItakuraM found, as in the XY 
case mentioned above, the absence of chiral fixed points 
of renormalization-group flow. However, he could not 
find numerical evidence of the first-order transition. He 
predicted that if the transition is of first order for the 
Heisenberg spins, it should occur at much larger lattice 
sizes which he was not able to perform at that time. En- 
couraged by the high performance of the Wang-Landau 
method, we decided to study the Heisenberg case in this 
work using the full Hamiltonian with very large lattice 
sizes. As shown below, we find indeed a first-order tran- 
sition in this case. 

The paper is organized as follows. Section II is devoted 
to the description of the model and the technical details 
of the Wang-Landau (WL) methods as applied in the 
present paper. Section III shows our results. Concluding- 
remarks are given in section IV. 

II. MONTE CARLO SIMULATION: 
WANG-LANDAU ALGORITHM 

We consider the stacking of triangular lattices in the z 
direction. The spins are the classical Heisenberg model 
of magnitude S = 1. The Hamiltonian is given by 

W = J^S,.S i , (1) 

where Si is the Heisenberg spin at the lattice site i, ^Zu j\ 
indicates the sum over the NN spin pairs Si and Sj both 



in the xy planes and in adjacent planes in the z direction. 
For simplicity, we suppose the same antiferromagnetic 
interaction J (J > 0) for both in-plane NN pairs and 
inter-plane NN ones. 

Recently, Wang and Landau 18 proposed a Monte Carlo 
algorithm for classical statistical models. The algorithm 
uses a random walk in energy space in order to obtained 
an accurate estimate for the density of states g(E) which 
is defined as the number of spin configurations for any 
given E. This method is based on the fact that a flat 
energy histogram H(E) is produced if the probability for 
the transition to a state of energy E is proportional to 
g(E)^ 1 . At the beginning of the simulation, the density 
of states (DOS) is set equal to one for all possible ener- 
gies, g(E) = 1. We begin a random walk in energy space 
(E) by choosing a site randomly and flipping its spin with 
a probability proportional to the inverse of the momen- 
tary density of states. In general, if E and E' are the 
energies before and after a spin is flipped, the transition 
probability from E to E' is 

p(E - E') = mm [g(E)/g(E'), 1] . (2) 

Each time an energy level E is visited, the DOS is mod- 
ified by a modification factor / > whether the spin 
flipped or not, i.e. g(E) — > g(E)f. At the beginning 
of the random walk, the modification factor / can be as 
large as e 1 ~ 2.7182818. A histogram H(E) records how 
often a state of energy E is visited. Each time the en- 
ergy histogram satisfies a certain "flatness" criterion, / 
is reduced according to / — * \ff and H(E) is reset to 
zero for all energies. The reduction process of the mod- 
ification factor / is repeated several times until a final 
value /fi na i which close enough to one. The histogram is 
considered as flat if 

H{E) > x%.(H(E)) (3) 

for all energies, where x% is chosen between 70% and 
95% and (H(E)) is the average histogram. 

The thermodynamic quantitie a 18 ! 19 can be evalu- 
ated by (E n ) = ±J2 E E n g(E) eM~E/k B T), C v = 

{E2 ^ } \ (M n ) = ±J2 E M n g(E)exp(-E/k B T), and 
X = { M 1~t*^ ) where Z is the partition function de- 
fined by I = J2E9( E ) ex P(- E / k BT)- The canonical 
distribution at any temperature can be calculated sim- 
ply by P(E,T) = ±g(E) exp(~E/k B T). 

In this work, we consider a energy range of interes t 20 ! 21 
(E m i n , E max ). We divide this energy range to R subin- 
tervals, the minimum energy of each subinterval is E l min 
for i — 1,2, ...,R, and maximum of the subinterval i is 
£Jj ax = E ^t , + 2AE, where AE can be chosen large 
enough for a smooth boundary between two subinter- 
vals. The Wang-Landau algorithm is used to calculate 
the relative DOS of each subinterval (E^^, E^^) with 
the modification factor /g na i = exp(10 -9 ) and flatness 
criterion x% = 95%. We reject the suggested spin flip 
and do not update g(E) and the energy histogram H(E) 
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of the current energy level E if the spin-flip trial would 
result in an energy outside the energy segment. The DOS 
of the whole range is obtained by joining the DOS of each 
subinterval (E l min + AE, £^ ax - AE). 



III. RESULTS 

We used the system size of N x N x N where N = 
72, 84, 90, 96, 108, 120 and 150. Periodic boundary 
conditions are used in the three directions. J = 1 is 
taken as the unit of energy in the following. 

The energy histograms for three representative sizes 
N = 96, N = 120 and N = 150 shown in Figs. CD H 
and [31 respectively. As seen, for N = 96, the peak is 
very broad, a signature of the beginning to of a double- 
maximum structure. The double peak begins really at 
TV = 120. We note that the distance between the two 
peaks, i. e. the latent heat, increases with increasing 
size and reaches 0.0025 for TV = 150. This is to be com- 
pared with the value ~ 0.009 for N = 120 in the XY 
ease———— Such a small value of the latent heat in 
the Heisenberg case explains why the first-order charac- 
ter was so difficult to be observed. For increasing sizes, 
the minimum between the peaks will be deepened to sep- 
arate completely the two peaks. Note that the double- 
peak structure is a sufficient condition, not a necessary 
condition, for a first-order transition. We give here the 
values of T c for a few sizes: T c = 0.95774, 0.95768 and 
0.957242 for A=96, 120 and 150, respectively. 

To explain why standard MC methods without his- 
togram monitoring (see for example Ref. [!) fail to see 
the first order character, let us show in Fig. 0]the energy 
vs T obtained by averaging over states obtained by the 
WL method for N = 96, 120 and 150. We see here that 
while the energy histograms show already a signature of 
double-peak structure at these big sizes, the average en- 
ergy calculated by using these WL histograms does not 
show a discontinuity: the averaging over all states erases 
away the bimodal distribution seen in the energy his- 
togram at the transition temperature. Therefore, care 
should be taken to avoid such problems due to averaging 
in MC simulations when studying weak first-order tran- 
sitions. 

Figures O and [HI show the magnetization and the sus- 
ceptibility for three sizes N = 96, 120 and 150. Again 
here, one does not see with one's eye the discontinuity 
of the magnetization at the transition even for N = 150. 
The averaging procedure erases, as for the energy, the 
detailed structure at the transition. 

At this stage it is interesting to make another check 
of the first-order character: in a first-order transition, 
the maximum of the susceptibility should scale with 
the system volume, namely N d where d is the system 
dimension.— We plot in Fig. [7]x ma;]: versus A in a In — In 
scale. The slope of the straight line is ~ 3.1 which is noth- 
ing but d within errors. This is a very strong signature 
of a first-order transition. 
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FIG. 1: Energy histogram for N = 96 at T c indicated on the 
figure. 
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FIG. 2: Energy histogram for N = 120 at T c indicated on the 
figure. 
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FIG. 3: Energy histogram for N = 150 at T c indicated on the 
figure. 
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FIG. 4: Energy versus T for N = 96, 120, 150. 
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becomes clearly of first-order only at a very large lattice 
size confirming the result of a nonperturbative RG cal- 
culations using an effective average Hamiltonian. The 
present work hence puts definitely an end to the long- 
standing controversial subject on the nature of the phase 
transition in Heisenbcrg STA. To conclude, let us empha- 
size that for complicated systems like this one, methods 



FIG. 5: Magnetization versus T for N =96, 120, 150. 




FIG. 6: Susceptibility versus T for N =96, 120, 150. 



IV. CONCLUDING REMARKS 

We have studied in this paper the phase transition in 
the Heisenberg STA by using the flat histogram tech- 
nique invented by Wang and Landau. The method is 
very efficient because it helps to overcome extremely long 
transition time between energy valleys in systems with a 
first-order phase transition. We found that the transition 
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FIG. 7: Maximum of susceptibility versus iV=96, 108, 120 
and 150 in the In — In scale. The slope is 3.1. See text for 
comments. 

well established for simple systems such as ferromagnets 
may encounter difficulties in dealing with the nature of 
the phase transition. Such difficulties can be solved only 
with high-performance MC simulations as the one used 
here, and a detailed analysis of the flow behavior as sug- 
gested by a nonperturbative RG calculation. 
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